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ABSTRACT 


We develop numerical procedures for constructing asymptotic solutions of 
certain nonlinear singularly perturbed vector two-point boundary value 
problems having boundary layers at one or both endpoints. The asymptotic 
approximations are generated numerically and can either be used as is or to 
furnish a general purpose two-point boundary value code with an initial 
approximation and the nonuniform computational mesh needed for such problems . 
The procedures are applied to a model problem that has multiple solutions and 
to problems describing the deformation of a thin nonlinear elastic beam that 
is resting on an elastic foundation. 
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Introduction 


Inxtxal value problems for stiff systems of ordxnary differentxal 
equatxons are now consxdered to be relatxvely tractxble numerically (cf. 
Enrxght et al. [7]). However, codes for stiff (or singularly perturbed) 
boundary value problems are not readily available, even though these problems 
arise in a great many applications. 

In this paper we consider asymptotic and numerical methods for singularly 
perturbed two-poxnt boundary value problems of the form 


X = f (x,y, t, e) 


ey = g{x,y, t, e) 


rs/ 


rs/ 


a(x(0) ,y(0) , e) = 0 , b(x( 1 ) ,y( 1 ) , e) = 0 , 


( 1 . 1 ) 


(1 ,2a,b) 


where x, y, a, and b are vectors of dimension m, n, q, and r = m + n - q, 

rw rsj /s/ 

respectively, and e is a small positive parameter. 

Although many special problems of this form can be solved by known 
asymptotic or numerical techniques, the general problem is very difficult and 
beyond our current understanding. The form of equations (1.1, 2) imply that, 
whenever g is not small, y varies rapidly relative to x. The behavior of the 
solution in these zones of rapid transition can be very complicated. For 
example, y can "jump" abruptly in a narrow boundary layer near t = 0 and/or 1 . 
These jumps can also occur at interior locations where solutions or their 
derivatives will become unbounded as e 0. The locations of the interior 
layers are generally unknown and must be determined as part of the solution 
process. Examples of these and other phenomena are discussed in, e.g., 
O'Malley [20], Kevorkian and Cole [18], Pearson [24, 25], Hemker [15], and 
Flaherty and O'Malley [11]. 
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The traditional numerical techniques for two-point boundary value 
problems all have difficulties with singularly perturbed problems unless the 
grid that is used for the discretization is appropriately fine, at least 
within boundary or interior layers. If the grid is not fine enough to resolve 
the layers, the computed solution typically exhibits spurious mesh 
oscillations. There are, however, special purpose schemes that can solve 
some singularly perturbed boundary value problems without using a fine 
discretization in transition regions. Most notable among these are the 
"upwind" or one-sided finite difference schemes (cf., e.g., Kreiss and Kreiss 
[19] or Osher [23]) and the exponentially weighted finite difference and 
finite element schemes (cf., e.g., Flaherty and Mathon [9] or Hemker[15]). 
These schemes must usually be either restricted to relatively simple problems 
or employ complicated algebraic transformations. 

In view of these theoretical and computational difficulties, we simplify 
problem (1.1, 2) considerably by assuming, in addition to natural smoothness 
hypotheses, that (i) g, a, and b are linear functions of the "fast" variable 
y, (ii) the n X n Jacobian 

G(x,t) ;= gy(x,y,t,0) (1.3) 

has a strict hyperbolic splitting with k > 0 stable and n - k > 0 unstable 
eigenvalues for all x and 0 < t <1, and (iii) q > k and r > n - k. A 
corresponding theory for problems with quadratic dependence on y is very 
limited (cf., e.g., Howes [17] which discusses second-order scalar equations). 
This, of course, limits extension of a numerical theory, but encourages 
further numerical experimentation. 

The assumed hyperbolic splitting restricts any rapid variations in y 
to occur in boundary layer regions near t = 0 and/or 1 . Thus, we 
unfortunately have eliminated many important and challenging problems having 
interior or "shock" layers. Some numerical work on these problems was done 
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by Kreiss and Kreiss [19], Osher [23], and O'Malley [22], 

In a series of three papers, Ascher and Weiss [2, 3, 4] show that 

symmetric, or centered, collocation schemes could be used on problems that 

satisfied assumptions similar to ours provided that appropriately fine meshes 

were used in the endpoint boundary layers. Our approach is somewhat different 

in that we use the assumed hyperbolic splitting to find an asymptotic solution 

of problem (1.1, 2) which is composed of a limiting outer solution (X (t), 

~0 

y (t)) and boundary layer corrections near t = 0 and 1. The limiting solution 
satisfies a reduced system, which is obtained from (1.1) by formally setting 
e to zero, i.e. , 

X = f(X ,Y ,t,0) , 0 = g(X ,Y ,t,0) . (1.4a,b) 

Because G is everywhere nonsingular, we can solve Eq. (1.4b) for 

rw 

Y = Y (X ,t) in a locally unique way, and there remains the m th order 
differential system (1.4a) for determining X (t). 


In order to completely specify the limiting solution, we must prescribe 
ra boundary conditions for Eq. (1.4a). We do this in Section 2 by providing a 
"cancellation law" that selects a combination of q - k initial conditions 
(1.2a) and of r - n + k terminal conditions (1.2b) to be satisfied by X (t). 
For more nonlinear problems, we note that such a cancellation law is much 
more difficult to specify (cf. O'Malley [21]). Boundary layer corrections are 
generally needed to compensate for the cancelled initial and ternanal 
conditions, and these are easily determined once X (t) and Y (t) have 


been found (cf. Section 2). 




In Section 3 we discuss a numerical procedure for calculating the 
asymptotic solution of Section 2. We implement the cancellation law by 
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using orthogonal transformations to reduce G(x(t),t) to a block triangular 
form with its stable and unstable eigenspaces separated. We also use the 
general purpose two-point boundary value code COLSYS to solve the reduced 
problem and then add numerical approximations of the boundary layer 
corrections. This approximation is considerably less expensive to obtain than 
solving the full stiff problem numerically and it has the advantage of 
improving in accuracy, without any additional computational cost, as the 
small parameter e tends to zero. However, when E is only moderately small, 
our asymptotic approximation may not be sufficiently accurate for some 
applications, so we have developed a procedure for generating an improved 
solution by using COLSYS to solve the complete problem (1.1, 2) with our 
asymptotic approximation as an initial guess. In order for this approach to 
succeed, we must also provide COLSYS with an initial nonuniform mesh that is 
appropriately graded in the boundary layers, and we give an algorithm for 
constructing such a mesh in Section 3. 

In Section 4 we apply our procedures to a third order model problem that 
has multiple solutions and to problems involving the deformation of a thin 
nonlinear elastic beam. These examples show that our methods can calculate 
accurate solutions of stiff problems for a very modest computational effort. 
While our algorithm for furnishing COLSYS with an initial guess and a 
nonuniform mesh does not seem to be optimal, it does offer some advantages 
over the more standard approach of continuation in e, where one starts with a 
large value of e (e.g., e = 1) and a crude initial guess of the solution and 
reduces e in steps so that the mesh is gradually concentrated into the 
boundary layer regions. 
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We also present two examples xn Sectxon 4 that are beyond the 
capabilities of our current methods because their solutions become unbounded 
as e ->-0. We include numerical results for these problems in this paper in 
order to show some of the several challenging effects that can occur with 
singularly perturbed problems. Finally, in Section 5, we discuss our results 
and present some suggestions for future investigations. 


2 . Asymptotic Approximation 

With the assumed hyperbolic splitting, we expect solutions of (1.1,2) to 
feature boundary layers in the fast y variable near both endpoints as e ->■ 0. 
Thus, It IS natural (cf. O'Malley [21]) to seek bounded uniform asymptotic 
expansions of the form 


x(t, e) =x(t,e) + eC( T, e) + enCa, e) "j 


y(t, e) = y(t, e) + p( t, e) + v( a, e) 


, 0 < t < 1 , 


( 2 . 1 ) 


where the outer solution (X(t,e), Y(t, e)) represents the solution 
asymptotically within (0,1), the initial layer correction (e5(T,e), u(t,e)) 
decays to zero as the stretched variable 


T = t/e 


(2.2a) 
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tends to xnfinity, and the terminal layer correction (en(cr/e), v(a.e)) 
goes to zero as the stretched variable 

a = (1 - t)/e (2.2b) 

approaches infinity. The outer solution and the boundary layer corrections 
are represented by expansions of the form 


- 



x(t, e) 


X (t) 
~3 

Y(t, e) 


Y (t) 
~3 

?( T, e) 

00 

- V 

K (T) 

p{ T, e) 

/N/ 

L 

3=0 

P (T) 

r) 

n( a, e) 


n ( o) 

v( a, e) 


V (a) 

n 


The limiting uniform approximation is obtained from (2.1) by letting e tend to 
zero, i.e.. 


x(t, e) = X (t) + 0(e), y(t, e) = Y (t) + p (t) + v (a) + 0(e) . (2.4) 

At t = 0 the fast vector y usually has a discontinuous limit, jumping from 

y(0,0) = Y (0) + p (0) to Y (0) at t = O'*’, An analogous Heaviside 
~ ~0 ~0 

discontinuity generally occurs near t = 1 . 
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The outer expansion (2.3a,b) must satisfy the full problem (1.1) within 

(0,1) as a power series in e; thus, the limiting solution (X ,Y ) will 

~0 ~0 

satisfy the nonlinear and non-stiff reduced system (1.4). As previously 

noted, since G(X ,t) (cf. Eq. (1.3)) is nonsingular we can solve Eq. (1.4b) 

~ HD 

for Y = Y (X ,t) in a locally unique way, so there remains the m th order 
hD HD HD 


nonlinear system (1.4a) for X (t). Later terms of the expansion (2.3a,b) 

HD 

satisfy linearized versions of the reduced system. For example, the 
coefficients of order e give 


X = fx(X ,Y ,t,0)x + f„(X ,Y ,t,0)Y + f (X ,Y ,t,0) , 

H Hs- HD H) ~1 H) HD ~1 ~e HD H) 

(2.4) 


Y = gx(X ,Y ,t,0)X + G(X ,t)Y + g (X ,Y ,t,0) . 

H) HD HD ~1 ~ HD ~1 ~£ HD HD 


We can determine Y (t) in terms of X , Y , and X from (2.4b) and, once again, 
~1 HD HD ~1 


there remains the m th order linear system (2.4a) for X . Similarly, for each 

~1 

j > 1 , we obtain a system of the form 


X = fx(X ,Y ,t,0)x + f„(X ,Y ,t,0)Y + a (X ,***,X ,t) , 

-'3 HD HD ^ H, HD HD ~j ~g-1 HD ~j-1 


(2.5) 

Y = gx(X ,Y ,t,0)X + G(X ,t)Y +3 (X ,***,X ,t) 

^_1 H) HD ~g ~ HD ~g ~3~1 HD ~j-1 


with successively determined inhomogeneous terms 
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In order to completely specify the outer expansion (2.3a,b), we must 

prescribe boundary conditions for the m-vectors X (t). Most critically, we 

need to specify m boundary conditions for the limiting slow vector X (t). 

M3 

It is natural to attempt to determine them by somehow selecting a subset of m 
combinations of the m + n boundary conditions (1.2) evaluated at e = 0. For 
scalar higher order linear differential equations, the first such 
"cancellation law" was obtained by Wasow [29] . Harris [14] obtained a more 
complicated cancellation law for linear systems with coupled boundary 
conditions and Ferguson [8] developed a numerical procedure for corresponding 
linear problems. These early works suggest that we should seek a cancellation 
law that ignores an appropriate combination of k initial conditions and of 
n-k terminal conditions. To this end, we must examine the boundary layer 
corrections and we begin by considering the initial layer correction (e5,p). 
Near t = 0, the terminal layer correction ( eg, v) may be neglected, so the 
representation of our asymptotic solution (2.1) requires the initial layer 
correction ( e?, u) to satisfy the nonlinear system 

d f/d T = dx/dt - dX/dt = f (x+ e^fY+y, ex, e) - f(X,Y,et,e) , 

rw ^ 

( 2 . 6 ) 

d y/d T = e(dy/dt - dY/dt) = g(x+e?,Y+y, et, s) - g(X,Y,ET,e) , 

on T > 0 and to decay to zero as t Substitution of the asymptotic 

expansion (2.3c, d) into (2.6) provides successive differential equations for 

the coefficients (5 ,y ). 

~3 

initial layer system 


In particular, when e = 0, we have the limiting 
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dC /dT = f(X (0),Y (0) + U ,0,0) - f(X (0),Y (0),0,0) , 

K3 ~ '-O 'K) ~ HD HD 

(2.7) 


dp /dT = g(X (0),Y (0) + p ,0,0) - g(X (0),Y (0),0,0) 


The decay requirement determines 


? ( t) = - / (d? (s)/dx)ds (2.8a) 

HD ~0 


as a functional of p , while p satisfies the conditionally stable system 

HD HD 


dp /dx = G(X (0),0)p . (2.8b) 

HD ~ HD HD 

We used (1.3) and the assumed linearity of g in y when obtaining (2.8b). 

/N# /N/ 

If g(x,y, t, e) were not linear in y, the initial layer correction would 

rw ^ 

satisfy a nonlinear differential equation which would generally be difficult 

to solve (cf. O'Malley [21]). Indeed, it would then be extremely difficult 

to specify what set of initial vectors p (0) would lead to decaying solutions 

HD 

of the boundary layer system (2.7b). Here, Eq. (2.4) is readily integrated to 
give 


G(X (0),0)x 

p ( x) = e 0 p (0) . (2.9) 

HD ~ 0 
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Thus, p will decay to zero as t -*■ " provided that 
■S3 


p (0) = P(X (0),0)p (0) , (2.10) 

'S3 ~ 'S3 ~0 


where P is a projection onto the k dimensional stable eigenspace of 

G(X (0),0). 

~ <S) 


Substituting (2.10) into (1.2a) and letting e -»■ 0, we see that the q 


limiting initial conditions take the form 


a(X (0),Y (0) + P(X (0),0)p (0),0) = 0 . (2.11) 

~ 'S3 'S3 ~ 'S3 »S3 ~ 


Now, using the linearity of a in y, we let 


A(x,t) = ^(x,y,t,0) (2.12) 


and further assume that A(X (0),0)P(X (0),0) has its full and maximal rank k. 

~ r-o -S) 


Then we can uniquely determine p (0) as a function of X (0) from k of the 

'S3 <S3 


equations (2.11). Having done this, initial conditions for the reduced 


problem can be determined from the remaining q - k conditions in (2.11). For 


the moment, we write these in the form 


$(X (0)) = 0 . (2.13) 

'S3 


In Section 3, we discuss a numerical procedure for determining P,p (0) and 

~ ~0 


$(X (0)). 
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The termxnal layer correction can be analyzed in an analogous manner. 


In particular, the leading term v ( a) satisfies 

~0 


G(X (1),1)a 

V ( a) = e~ V (0) . (2.14) 

~0 ~0 


Now, V will decay to zero as a ® provided that 
-0 


V (0) = Q(X (1),1)v (0) , (2.15) 

~0 ~ K ) ~0 


where Q is a projection onto the n - k dimensional unstable eigenspace of 


G(X (1),1). Substituting (2.15) into (1.2b) and letting e -*■ 0 gives 
~ ^>0 


the r limiting terminal conditions as 


b(x (1),Y (1) + Q(X (1),1)V (0),0) = 0 . (2.16) 


We let 


B(x,t) = by(x,y,t,0) (2.17) 

^ rsr /V /s/ 


and assume that B(X (1),1)Q(X (1),1) has full rank n - k. Then we can solve 
~ ~0 ~ ^ 


(2.17) for V (0) and the remaining r - n + k conditions specify terminal 
'K) 

conditions for the limiting problem, which we denote by 


'?(X (D) = 0 . (2.18) 

0 


The reduced problem consists of the nonlinear reduced differential 


equation and the m separated nonlinear boundary conditions (2.13, 18). If 


it IS solvable, it may have many solutions; however, corresponding to any of 


Its isolated solutions (X (t),Y (t)), one can expect to find a solution of 

the original problem (1.1, 2) that converges to (X (t),Y (t)) on 0 < t < 1 

~0 M) 
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as e ^ 0. Sufficient hypotheses to obtain an asymptotic solution having the 
form of (2.1) are provided by Hoppensteadt [16] and others. For this reason, 
we shall merely indicate the considerations that are involved in obtaining 
further terms in the initial and terminal layer expansions and boundary 
conditions for the outer expansion. 

Additional terms of the initial layer expansion (2.3c,d) are determined 
by equating the coefficients of in the nonlinear system (2.7), i.e. 


d? /dx = f„(X (0),Y (0) + u (t),0,0)p + Y (t) , 

~0 ~0 ~0 '^-l 

(2.19) 


d p /d X = G( X ( 0) ,0) y +5 ( x) , 

~ ~0 'YJ '^-1 


for 3 > 1 , where the inhomogeneous terms are exponentially decaying as x + “ 

because ? and p , i = 1, ..., 3 -I , and their derivatives behave in 

'^-1 -XL-1 

this manner. The linear system (2.19) may Jdo integrated to yield 


5 ( x) = -/ (d5 (s)/dx)ds , 
n X 

( 2 . 20 ) 

G(X (0),0)x X G(X (0),0)(x-s) 

p ( x) = e~ p (0) + / e ® 6 (s)ds 

no 


We see that 5 ( x) decays as x increases and p (x) will decay provided that 

n n 


p (0) lies in the unstable eigenspace of G(X (0),0), i.e., 
~] ~ ~0 


p (0) = P(X (0),0)p (0) 
~g ~ ~0 ~3 


( 2 . 21 ) 
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Using (2.1) and (2.3a,b) we find that the coefficient of e in the 


initial condition ( 1 . 2 a) has the form 


a (X (0),Y (0) + y (0),0)X (0) + A(X (0),0)[Y (0) 

K ) ~~0 

( 2 . 22 ) 


+ P(X (0),0)y (0)] = C 
~ ”^-1 


Since A(X (0),0)P(X (0),0) has its maximal rank k, we can determine y (0) 
~ ~ 'K) ~j 


from k of these equations, and the remaining q - k equations determine linear 


equations for X (0). The situation for the terminal layer correction is 
completely analogous; thus, v (0) and the terminal conditions for X (1) are 

n ~j 


determined from linear equations of the form 


bjj(X (1),Y (1) + V (0),0)X (1) + B(X (1),1)[Y (1) + Q(X (1),1)v (0)] = 0 

'O ~~0 ~~0 ~3 ~]-1 


To summarize, we have shown that the j th (3 >1) term in the outer 

expansion satisfies cin m th order linear boundary value problem consisting of 

Eq. (2.5) and a set of m linear boundary conditions determined from (2.22) 

and (2.23). It is a linearization of the problem for X (t). 

~0 


(2.23) 
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3 . Numerical Procedure . 

In this section we discuss a numerical procedure for finding the limiting 
uniform asymptotic solution (2.4). It consists of solving the limiting outer 
problem (1.4, 2.13, 2.18) and determining boundary layer corrections from 
(2.9) and (2.14). 

Our first task is to find the projections P and Q and we do this by 
finding the Schur decomposition of the matrix G at t = 0 and t = 1 . In 
particular, at t = 0 we find an orthogonal matrix E(x(0),0) such that 


T (x(0),0) 


G(x(0) ,0)E(x(0),0) = E(x(0),0) 


0 


U(x(0),0) 


T (x(0),0) 


(3.1 ) 


where T is k x k and upper triangular with the stable eigenvalues of 

G(x(0),0), and T is upper triangular with the remaining n - k unstable 
~ ~ 

eigenvalues. The decomposition (3.1) can often be obtained analytically; 
however, when this is not possible or practical it can be determined 
numerically by using the QR algorithm (cf. Golub and Wilkinson [13], 

Ruhe [26], and Bjork [5] for specific procedures). 

We partition E after its k th column as 

E = [E "e ] (3.2) 


and note that E spans the stable eigenspace of G at t = 0 and 


T 

P = E E 


(3.3) 


IS the desired projection onto this eigenspace 
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Substituting (3.3) into (2.11) gives 


T 

a(X (0), Y (0) + E (X (0),0)E (X (0),0)y (0),0) = 0 . 


as the equation for determining u (0) and $(X (0)). Since 

~0 ~ 


A(X (0),0)E (X (0),0) is of rank k, we construct a q x q matrix 
~ 'O -MD 


T T _T 

L = [L L ] 


(3.4) 


(3.5a) 


that reduces it to echelon form, i.e.. 


- 



L 

A(X (0),0)E (X (0),0) = 

V 

L 


0 





(3.5b) 


where V is k x k and nonsingular. Multiplying Eq. (3.4) by L, using the 


linearity of a in y, and Eq. (3.5) gives the initial layer jump p (0) and the 
~ ~ ~0 


q - k initial conditions (2.13) for the reduced problem, respectively, as 


-1 

u (0) = -E (X (0),0)V L a(X (0),Y (0),0) , 

(3.6) 


$(X (0)) := L a(X (0),Y (0),0) = 0 . 


We find the terminal layer jump and the r - (n-k) terminal conditions 


for the reduced problem in an analogous fashion with the exception that we 


define E(x(1),1) such that 
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G(x(1),1 )E(x(1),1) = E(x(1),1) 


T (x(1),1) 
0 


U(x(1),1) 


T (x(1),1) 


(3.7) 


which we partition after its (n - k) th column as 


E = [E E ] 

»m| 


(3.8) 


In parallel with Eqs. (3.1) and (3.2), the matrices T , T , and E contain 

rww '*'■*4* 

the )<: stable eigenvalues, the n - )c unstable eigenvalues, and span the 
unstable eigenspace, respectively, of G at t = 1 . Our reasons for switching 
the positions of the matrices containing the stable and unstable eigenvalues 
of G IS that we are unaware of a simple and stable computational procedure 
for finding a set of vectors that span a given subspace and are not in the 
leading columns of an orthogonal matrix like E (cf. Golub and Wilkinson [13]). 
Now, following the procedure that we used for the initial layer, we take 

T 

Q(X (1),1) = E (X (1),1)E (X (1),1) (3.9) 

as our projection onto the (n - k) dimensional unstable eigenspace of 

G(X (1),1) and construct an r x r matrix 
~ -0 

T T _T 

R = [R R ] (3.10a) 


that reduces the rank n - k matrix B(X (1),1)E (X (1),1) to echelon form, 

~ ~0 


r n 


m ” 

R 


V 




R 

B(X (1),1)E (X (1),1) = 
~ K) M3 

0 



- ~ - 


(3.10b) 
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where V is (n-k) x (n-k) and nonsingular. Multiplying Eq. (2.16) by R 
and using Eqs. (3.9) and (3.10), we find the terminal layer jump and terminal 
conditions for the reduced problem, respectively, as 

-1 

V (0) = -E (X (1),1)V R b(X (1),Y (1),0) , 

“SD ~0 "H- M3 

(3.11 ) 

'KX (D) := R b(X (1),Y (1),0) = 0 
^ ^ /»-0 ^ 


Since the reduced problem (1.4), (3.6b), and (3.11b) is not stiff, we can 

use any good code for two-point boundary value problems (cf. Childs et al.[6]) 

to solve It, and we have chosen to use the collocation code COLSYS of Ascher, 

Christiansen, and Russell [1]. The reduced problem is generally nonlinear and 

since COLSYS solves nonlinear problems using a damped Newton method, we have 

to supply formulas for evaluating the Jacobians of f, Y, 3>, and f with 

^ 

respect to X. We do this, but introduce an error, by providing analytical 
formulas for these Jacobians that neglect the influence of the derivatives of 
E, L, R, and G. (These derivatives will be small when the related 

/>w 

subspaces are nearly constant). This procedure failed to converge once on 
Example 1 of Section 4 and a minor modification to the Jacobian of $ restored 
convergence; however, an alternative possibility would be to approximate the 
Jacobians by finite differences. 

We start the Newton iteration with a uniform mesh and an initial guess 

( 0 ) 

X (t) for X (t). In section 4, we used the default initial guess that is 
M) M) 

provided by COLSYS for Example 2 and a constant initial guess for Example 1 . 
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This latter choice was necessary as Example 1 has three solutions. At each 

(P) 

Iteration step, we calculate an approximation E(X (t),t) to E(X(t),t) 

0 

(P) 

for t = 0 and 1 as the Schur decomposition of G(X (t),t). The examples 

of Section 4 were calculated using analytical formulas for E rather than the 

numerical procedures of Golub and Wilkinson [13], Ruhe [26], or Bjork [5]. 
(p) (p) 

Finally, L and R are obtained by using Gaussian Elimination to row 


(p) (p) (p) (p) 

reduce A(X (0),0)E (X (0),0) and B(X (1),1)E (X (1),1), 

respectively. 

When this procedure converges to (X (t),Y (t)), we calculate boundary layer 

~0 

corrections u ( t) and u (a),, for a given value of e, using Eqs. (2,9). (3.6a), 
(2.14), and (3.11a), and add these to the reduced solution in order to get 
the 0(e) asymptotic approximation (1.4), For moderately small values of e, 
this approximation may not provide a sufficiently accurate representation of 
the solution and, in this case, we use it as an initial guess to COLSYS and 
solve the complete problem (1.1, 2). However, this procedure may fail 
unless we also provide COLSYS with an initial nonuniform partition 


ir:={0=t <t <***<t = l} 

0 1 N 


(3.12) 


that IS appropriately graded within the boundary layers. Following Ascher, 
Christiansen, and Russell [1], we seek to find ir such that the error on each 
subinterval satisfies 


< 5(1 + 



1 


1 , 2 , ..., 


N 


(3.13) 
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where 6 is a prescribed tolerance. 


T T T 
u = [x , y ] , 


(3.14) 


e(t) is the difference between u(t) and its collocation approximation. 


|u| I = max |u(t)| , and |u(t) 
~ i t <t<t ~ 
i -1 i 


max |u (t) 
1 <3 <m+n 3 


(3.15) 


We assume that the final partition selected by COLSYS to solve the 
reduced problem satisfies (3.13) outside of boundary layer regions and 
we seek to refine it within the boundary layers. We further assume that de- 
rivatives of u can adequately be approximated by either p (x) or v (a) in the 
~ ~0 ~0 

left or right boundary layer, respectively. 

It is known (cf. Russell and Christiansen [27]) that if the solution of 
(1.1,2) IS smooth 


I II II ^3+1 ) I , 3+1 3+2 

|e|| =c|ju |jh + 0 (h ) 

~ 1 3 ~ 


(3.16) 


for collocation at the image of 3 Gauss-Legendre points per subinterval. 

Here c is a known constant, 

3 


h = t - t , and h = max h 
1 1 1-1 Ki<N ^ 


(3.17) 


In the left boundary layer we approximate u in (3.16) by p using (2.9) and 
attempt to find a partition that satisfies 


c h^ (t/E)|j « 6(1 + ||u|| ) . 

3 1 ~0 1 1 


(3.18) 
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Finally, we use (2.9) and (3.1) to approximate p and the subinterval lengths 


t - t 
1 1-1 


e 


(-) 

a 


5(1 + Null ) 
~ 1 

c I p (t /e) 1 

j ~0 1-1 


■ 1/(D + 1) 


(3.19) 


where cl is the magnitude of the largest diagonal element of T (X (0),0). 

0 

A similar formula can be obtained for selecting subintervals in the right 
boundary layer. 

Starting with i = 1 , we use Eq. (3.19) to generate a partition until we 
either reach t = 1/2 or a point where a subinterval length selected by Eq. 
(3.19) is larger than that used locally by COLSYS to solve the reduced 
problem. We then repeat the procedure in the right boundary layer. 

We have written a computer code called SPOOL that implements the 
algorithms that are described in this section; thus, it (i) uses COLSYS to 
solve the reduced problem, (ii) calculates and adds appropriate boundary layer 
corrections to the reduced problem, and (lii) (optionally) suggests a mesh 
that can be used by COLSYS to solve the complete problem. 


4 . Examples . 

In order to appraise the performance of SPCOL, we have applied it to a 
problem involving the deformation of a thin nonlinear elastic beam (Example 1) 
and a third order model problem that has multiple solutions (Example 2). 
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Example 1 . We consider problems involving the deformation of a nonlinear 


elastic beam that is resting on an elastic foundation with unit spring 
constant and is subjected to the combined action of a horizontal end thrust P 
and a unit uniform lateral load. This problem is discussed in detail in 
Flaherty and O'Malley [11] and herein we only present the governing 

differential equations, which in dimensionless form are 

• • • 

X = cos X , X = sin X , X = y , (4,1a,b,c) 

1 3 2 3 3 1 


• • 

ey = -y , ey = (x -1)cos x - Ty , (4.1d,e) 

1 2 2 2 3 1 


where 


T = sec X + ey tan x . (4. If) 

3 2 3 


The slow variables (x ,x ) and x represent the Cartesian coordinates and the 

1 2 3 

tangent angle of a material particle on the centerline of the beam that was at 


the Cartesian location (t,0) in the undeformed state. The fast variables y 

1 

and y are the internal bending moment and transverse shear force, 

2 

respectively. Finally, the small parameter is 

2 2 

E = EI/PL , (4.2) 


where El is the flexural rigidity and L is the length of the beam; thus, our 


beam is much stronger in extension than it is in bending 
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This example does not precisely fit our hypotheses since the axial force 

T IS a function of the fast variable y and, thus, g„ also depends on y. 

2 ~ 

However, our theory and methods will still apply as long as y remains bounded 

and lx I < Tr/2 as e 0. Flaherty and O'Malley [11] show that unbounded 
3 

solutions can occur when certain types of boundary conditions are prescribed 
for Eq. (4.1). In this paper we present results for the following three sets 
of boundary conditions: 

(i) . X (0, e) = X (0,e) = y (0,e) = x (1,e) = y (1,e) = 0 , (4.2a) 

12 12 1 

(ii) . X (0, e) = 0, -lOx (0,e) + y (0, e) = 0, -x (0,e)+10y (0,e) = 0 

■■ 2 2 3 1 (4.2b) 

lOx (1,e) + y (1,e) = 0, lOx (1,e)+y (1,e) = 0 , 

2 2 3 1 

(ill). X (0, e) = X (0,e) = X (0,e) = x (1,e) = x (1,e) = 0 . (4.2c) 

1 2 3 2 3 

Equations (4.2a) correspond to "simple supports", Eqs. (4.2c) correspond to 
"clamped supports", and Eqs. (4.2b) correspond to elastic supports that are 
almost simply supported at t = 0 and almost clamped at t = 1 . Conditions 
(4.2b) could arise because, say, friction introduces some coupling between 
lateral and rotational effects at the supports. As we shall see, y remains 
bounded for conditions (4.2a,b), but becomes unbounded as e ->■ 0 when 
conditions (4.2c) are applied. The problem is that the boundary conditions 
for the clamped beam only involve the slow variables and the slow vector x 
cannot generally satisfy all five of them without having boundary layers. 

This in turn forces the fast vector y to become unbounded like 0(l/e) at 
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the endpoints. Thus, the solution cannot have an asymptotic expansion of 

the form of Eq. (2.1); however, an appropriate asymptotic expansion was 

obtained by Flaherty and O'Malley [11]. We do not repeat those results here, 

but in order to emphasize the diverse behavior that can occur with nonlinear 

singularly perturbed problems, we present solutions for x , x , and y 

2 3 3 

corresponding to each of the boundary conditions (4.2a), (4.2b), and (4.2c) 
in Figures 1 , 2, and 3, respectively. 

Our methods apply to problems having boundary conditions (4.2a) and 
(4.2b) and, in these cases, the orthogonal matrix 


2 - 1/2 

E(x(0),0) = (1+a ) 



(4.3a) 


where 


2 

a = sec x (0) (4.3b) 

3 


reduces 


0 


G(x(0),0) = 


-1 



0 J 


(4.4) 


T 

to the Schur form given by (3.1) at t = 0 while E (x(1),1) will reduce 
G(x(1),1) to the form given by (3.7) at t = 1 . 

We solved Eq. (4.1) with conditions (4.2a) and (4.2b) in two ways: 


(i) using COLSYS to solve the complete problem with continuation from a large 



-1.00 -0.80 -0.60 -0.40 -0.20 -0.00 0.20 0.00 0.04 0.08 0.12 0.16 0.20 


2h. 



Figure 1. Nimerical solution of simply supported beam. Example 1 with 
boundary conditions given by equation (i+.2a). 








2.00 



Figure 2. Numerical solution of elastically supported beam. Example 1 
with boundary conditions given by equation (i+.2b) . 
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O 






Figure 3. Numerical solution of clamped beam. Example 1 with boundary- 
conditions given by equation (4.2c). Note that y and y 
are multiplied by e. 12 
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to a small value of e and (ii) using our code SPCOL to compute an initial 
asymptotic approximation and to recommend a nonuniform mesh and using this 
with COLSYS to calculate an improved solution. All calculations were 
performed using double precision arithmetic on an IBM 3033 computer, two 

collocation points per subinterval, and an error tolerance 6 (cf. Eq. (3.19)) 

-6 -3 

of 10 for slow variables and 10 for fast variables. 

Our results for the normalized CP times and the number of subintervals 

(NSUB) that are either used by COLSYS or recommended by SPCOL are shown in 

Tables 1 and 2 for the simply supported beam and in Tables 3 and 4 for the 

elastically supported beam. Tables 1 and 3 contain the continuation results 

and Tables 2 and 4 contain the SPCOL results with COLSYS improvement. The CP 

times (for all examples) were normalized with respect to the e sequence in 

Table 1 . Differences between our initial asymptotic approximation and the 

final solution obtained by COLSYS are shown for x (1/2, e) and y (0,e) for the 

2 2 

simply supported beam in Table 5 and for x (0,e) and y (0,e) for the 

3 2 

elastically supported beam in Table 6. All of the differences are decreasing 

like 0( e) as expected. Differences that are recorded as zero are less than 
-8 

10 . 

The results reported in these Tables need some additional explanation. 

The number of subintervals and CP times used with continuation depended 
heavily on the e sequence that was used. The results in Tables 1 and 3 are 
about the best insofar as they gave the smallest total CP time for the 
sequence. We see in almost every instance that the COLSYS correction is using 
about twice the number of subintervals that were suggested by SPCOL. This 
mesh doubling strategy is often used in COLSYS to estimate errors or when the 
Newton iteration has convergence difficulties. Thus, in some sense our mesh 
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£ 

NSUB 

CP 

TOTAL CP 

-1 

10 

80 

6.1 

6.1 

-2 

10 

72 

6.3 

12.5 

-4 

10 

112 

18.4 

30.9 

-6 

10 

158 

27.2 

58.1 

-8 

10 

254 

41 .9 

100.0 


TABLE 1 . EXAMPLE 1 WITH SIMPLE SUPPORTS . NUMBER OF SUB- 
INTERVALS (NSUB) AND CP TIMES TO SOLVE THE PROBLEM 
BY COLSYS WITH CONTINUATION IN e. TOTAL CP IS THE 
ACCUMULATED TIME FOR THE £ SEQUENCE. 


£ 

SPCOL 

COLSYS 

CORRECTION 



NSUB 

CP 

NSUB 

CP 

TOTAL CP 

-1 

10 

20 


1 .3 

80 

5.7 

7.0 

-2 

10 

28 


1 .3 

112 

8.7 

10.0 

-4 

10 

34 


1 .3 

136 

9.0 

10.3 

-8 

10 

35 


1 .3 

92 

9.3 

10.6 


TABLE 2. EXAMPLE 1 WITH SIMPLE SUPPORTS. NUMBER OF SUB- 
INTERVALS (NSUB) AND CP TIMES TO SOLVE THE PROBLEM BY 
SPOOL AND TO IMPROVE IT BY COLSYS. THE CP TIMES 
INCLUDE THE TIME TO CALCULATE THE REDUCED SOLUTION, 
WHICH WAS 1.3 TIME UNITS. TOTAL CP IS THE SUM OF THE 
SPCOL CP AND THE COLSYS CP. 














e 


NSUB 


CP 


TOTAL CP 


-1 




10 

80 

6.9 

6.9 

-2 




10 

78 

6.3 

14.6 

-4 




10 

78 

16.8 

31 .4 

-6 




10 

156 

38.3 

69.7 

-8 




10 

100 

16.4 

86.1 


TABLE 3. EXAMPLE 1 WITH ELASTIC SUPPORTS. NUMBER OF SUB- 
INTERVALS (NSUB) AND CP TIMES TO SOLVE THE PROBLEM 
BY COLSYS WITH CONTINUATION IN £. TOTAL CP IS THE 
ACCUMULATED TIME FOR THE e SEQUENCE. 



SPCOL 

COLSYS 

CORRECTION 


e 

NSUB 

CP 

NSUB 

CP 

TOTAL CP 

-1 

10 

40 

3.9 

100 

10.2 

14.1 

-2 

10 

47 

3.9 

94 

10.5 

14.4 

-4 

10 

56 

3.9 

112 

12.8 

16.7 

-8 

10 

57 

3.9 

134 

16.8 

20.7 


TABLE 4. EXAMPLE 1 WITH ELASTIC SUPPORTS. NUMBER OF SUB- 
INTERVALS (NSUB) AND CP TIMES TO SOLVE THE PROBLEM 
BY SPCOL AND TO IMPROVE IT BY COLSYS. THE CP TIMES 
INCLUDE THE TIME TO CALCULATE THE REDUCED SOLUTION, 
WHICH WAS 3.8 TIME UNITS. TOTAL CP IS THE SUM OF THE 
SPCOL CP AND THE COLSYS CP. 
















e 

Ax (1/2, e) 
2 

Ay (0,e) 
2 

-1 

-3 

-2 

10 

7.1x10 

3.2x10 

-2 

-5 

-3 

10 

6.7x10 

3.6x10 

-4 


-5 

10 

0 

3.6x10 

-8 



10 

0 

0 


TABLE 5. EXAMPLE 1 WITH SIMPLE SUPPORTS. DIFFERENCES BETWEEN 
SPOOL AND COLSYS SOLUTIONS, x.e., A( ) := | ( )sPCOL 
~ ( ) COLSYS I* 


e 

Ax (0,e) 
3 

Ay (0,e) 
2 

-1 

-1 

-2 

10 

1 .3x10 

4.2x10 

-2 

-2 

-3 

10 

1 .4x10 

5.2x10 

-4 

-4 

-5 

10 

1 .5x10 

5.4x10 

-8 



10 

0 

0 


TABLE 6. EXAMPLE 1 WITH ELASTIC SUPPORTS. DIFFERENCES BETWEEN 
SPOOL AND COLSYS SOLUTIONS, i.e., A( ) := |( ) SPOOL 
“ ( ) COLSYS I* 
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strategy is doing about as well as can be expected; however, it seems that 
fewer points should be necessary. We tried placing the subintervals according 
to a pointwise error criteria, as suggested by Ascher and Weiss [2, 3, 4], 
rather than the global criteria used in Eq. (3.19), but this gave very similar 
results (cf. Flaherty and O'Malley (12]). We also tried suggesting an initial 
mesh to COLSYS that consisted of every other point in the mesh suggested by 
SPOOL. This IS clearly a risky strategy, since collocation at the 
Gauss-Legendre points is known to produce oscillations unless the mesh is 

appropriately fine in the boundary layers (cf. Ascher and Weiss [2]). 

-8 

Nevertheless, this did give some improvement for values of E > 10 (cf. 
Flaherty and O'Malley [12]). Perhaps the results could be improved further by 
using higher order collocation and/or collocation at the Gauss-Lobatto points 
as suggested by Ascher and Weiss [2, 3, 4]. 

-8 

We see from Tables 1 to 4 that for e = 10 the SPOOL solution can be 
computed in less than 5% of the time of the continuation solution and the 
COLSYS improvement with the SPOOL solution as an initial guess can be computed 
in less than 24% of the time of the continuation solution for both simple and 
elastic supports. 


Example 2 . We consider the third order model problem 


X 


1 


- X f 




2 

= a (x)y 

1 


+ 8x(1-x) 


(4.5a,b,c) 


with 


a(x) = 1 + 2x 


(4.5d) 
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and the linear boundary conditions 


x(0, e) + y (0,e) = 0, -Yx(0,e) + y (0,£) = 0, x(1,e) + y (1,e) = 0. (4.6) 

1 2 1 


The matrix G(x,t) for this example is the negative of that given by (4,4) for 


Example 1 with a now being given by (4.5d). Thus, G has one negative and one 

positive eigenvalue provided that a(x) is nonzero and G may be reduced to 

T 

Schur form at t = 0 using the orthogonal matrix E (x(0),0) and at t = 1 using 
E(x(1),1) (with E(x,t) given by Eq. (4.3a)). 


Flaherty and O'Malley [10] studied this problem and showed that the 


reduced system is 


. 2 

X=1-X, Y =0, a(X)Y +8X (1-X ) = 0 (4.7) 

0 0 20 0 10 0 0 


with the initial condition 


1 o(X (0))|[X (0) + Y (0)] - YX (0) = 0 . (4.8) 

0 0 10 0 


They show that there are three solutions of (4.7,8) for each value of the 

constant y pi^ovided that there are no "turning points", i.e., provided that 

there are no values of x(t) for which a(x) = 0 on 0 < t < 1 . The three 

solutions can be characterized by their value of X (0) which is determined as 

0 

X (0) = 0 , -[ ys- 6±/( ys-4) ^ + 48] , s = sgn(a(X (0)) . 

0 4 0 


(4.9) 
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For y = 2 the three values of X (0) are 0, 0.803, and -4.29 and the three 

0 

corresponding solutions for y (t, e) are shown in Figure 4. For X (0) = 0, the 

1 0 

initial layer correction y ( x) is trivial; however, the other two solutions 

0 

have initial layer jumps. 

It can be easily verified that a(X (t)) has a zero on 0 < t < 1 when 

0 

(-3.08 «)-3e/2 +1 < X (0) < -1/2. In this case (4.5) has a turning point and 

0 

Y becomes unbounded. Our theory and methods do not apply in this case; 

10 

however, if e is not too small, the solution of (4.5) can be calculated using 

COLSYS. In order to contrast solutions with and without turning points, we 

illustrate y (t, e) for y = -2 and X (0) = -2.80 in Figure 5. 

1 0 

Solutions obtained using SPOOL and the corresponding COLSYS corrections 

are shown for y = 2 and X (0) = 0, 0.803, and -4.29 in Tables 7, 9, and 11, 

° -6 
respectively. The COLSYS correction failed to converge for e < 10 when 

X (0) =0 and -4.29. We have no explanation as to why the solution with 
0 

X (0) = 0.803 was so much easier to calculate. The relative difference 
0 

between the SPCOL and COLSYS solutions for x(l, e) and y (1 ,e) are shown in 

2 

Table 13 for T = 2 and X (0) = -4.29. These results are typical of those 

0 

that we obtained for all three solutions. 

Using COLSYS with continuation in E and the default initial guess can find 

at most one solution, and, for this example, it found the X (0) =0 solution. 

0 

The results of this calculation are shown in Table 8 for y = 2. Although sev- 

-6 

eral e sequences were tried, we were unable to obtain convergence for e < 10 
Again, this situation could possibly be remedied by using collocation at 
the Gauss-Lobatto points as in Ascher and Weiss (2, 3, 4]. The other two 
solutions when T = 2 can also be calculated using continuation in e provided 
that we use a suitable initial guess. Results for the solutions corresponding 
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to X (0) = 0.803 and -4.29 are presented in Tables 10 and 12, respectively, 

0 

using continuation with SPCOL furnishing an initial guess. These results seem 
to point to the possibility of using a combination of asymptotics and 
continuation to solve singular perturbation problems. 


5 . Discussion . 

We have obtained asymptotic approximations for a restricted class of 
nonlinear singularly perturbed two-point boundary value problems and have 
shown how to construct approximate solutions numerically and use them to 
suggest a nonuniform mesh that may be used as input to a two-point boundary 
value code in order to calculate improved solutions. Clearly this approach 
offers some advantages over the more standard technique of continuation in £ 
steps; however, the picture is far from clear and several questions still 
remain as to how best to use asymptotic analysis in conjunction with numerical 
analysis . 

In Example 2 of Section 4 we have shown that asymptotic methods may be 
used to distinguish different solutions in problems having multiple solutions. 
These asymptotic approximations may be used to provide initial guesses to a 
two-point boundary value code. 

In Example 1 of Section 4 we have shown that unbounded solutions can 
result from seemingly minor changes in the boundary conditions of singularly 
perturbed boundary value problems. Other very diverse behaviors can occur 
when turning point problems are considered (cf., e.g., Kevorkian and Cole [18] 
or O'Malley [20]). Since phenomena cannot easily be predicted, a sensible 
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SPOOL 

COLSYS 

CORRECTION 


G 

NSUB 

CP 

NSUB 

CP 

TOTAL CP 

-1 

10 

40 

0.6 

88 

6.3 

6.9 

-2 

10 

44 

0.6 

88 

6.2 

6.8 

-4 

10 

47 

0.6 

192 

18.2 

18.8 

-6 

10 

47 

0.6 


failed 



TABLE 7. EXAMPLE 2 WITH Y = 2 AND XgCO) =0. NUMBER OF SUBINTERVALS 
(NSUB) AND CP TIMES TO SOLVE THE PROBLEM BY SPOOL AND TO 
IMPROVE IT BY COLSYS. THE CP TIMES INCLUDE THE TIME TO 
CALCULATE THE REDUCED SOLUTION, WHICH WAS 0.5 TIME UNITS. 
TOTAL CP IS THE SUM OF THE SPOOL CP AND THE COLSYS CP. 


e 

NSUB 

CP 

TOTAL CP 

-1 

10 

40 

1 .8 

1 .8 

-2 

10 

44 

3.3 

5.2 

-4 

10 

264 

13.4 

18.6 

-5 

10 

372 

20.2 

38.7 

-6 

10 

___j 


failed 



TABLE 8. EXAMPLE 2 WITH Y = 2 and XqCO) = 0. NUMBER OF SUBINTERVALS 
(NSUB) AND CP TIMES TO SOLVE THE PROBLEM BY COLSYS WITH 
CONTINUATION IN e FROM e = 10"”' . THE DEFAULT INITIAL GUESS 
THAT IS PROVIDED IN COLSYS WAS USED TO START THE CONTINUATION' 
SEQUENCE. TOTAL CP IS THE ACCUMULATED TIME FOR THE 
SEQUENCE. 
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e 

SPOOL 

COLSYS 

CORRECTION 

TOTAL CP 

NSUB 

CP 

NSUB 

CP 

-1 






10 

42 

1 .5 

42 

3.0 

4.5 

-2 






10 

52 

1 .6 

52 

3.0 

4.6 

-4 






10 

57 

1 .6 

58 

2.6 

4.2 

-6 






10 

57 

1 .6 

114 

10.9 

12.5 


TABLE 9. EXAMPLE 2 WITH Y = 2 AND Xq = 0.803. NUMBER OP SUBINTERVALS 
(NSUB) AND CP TIMES TO SOLVE THE PROBLEM BY SPOOL AND TO 
IMPROVE IT BY COLSYS. THE CP TIMES INCLUDE THE TIME TO 
CALCULATE THE REDUCED SOLUTION, WHICH WAS 1.5 TIME UNITS. 
TOTAL CP IS THE SUM OF THE SPOOL CP AND THE COLSYS CP. 


e 

NSUB 

CP 

TOTAL CP 

-4 

10 

58 

2.6 

2.6 

-5 

10 

58 

2.4 

5.0 

-6 

10 

70 

4.3 

9.3 


TABLE 10. EXAMPLE 2 WITH Y = 2 AND Xq(0) = 0.803. NUMBER OF 

SUBINTSRVALS (NSUB) AND CP TIMES TO SOLVE THE PROBLEM 
BY COLSYS WITH CONTINUATION IN e FROM E = 10“'^. TOTAL 
CP IS THE ACCUMULATED TIME FOR THE SEQUENCE. 
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£ 

SPOOL 

COLSYS 

CORRECTION 

TOTAL CP 

NSUB 

CP 

NSUB 

CP 

-1 






10 

44 

0.9 

62 

3.6 

4.5 

-2 






10 

52 

0.9 

84 

3.8 

4.7 

-4 






10 

59 

0.9 

232 

15.3 

16.2 

-6 






10 

59 

0.9 


failed 



TABLE 11. EXAMPLE 2 WITH Y = 2 AND Xq(0) = -4.29. NUMBER OF 

SUBINTERVALS (NSUB) AND CP TIMES TO SOLVE THE PROBLEM 
BY SPOOL AND TO IMPROVE IT BY COLSYS. THE CP TIMES 
INCLUDE THE TIME TO CALCULATE THE REDUCED SOLUTION, WHICH 
WAS 0.9 TIME UNITS. TOTAL CP IS THE SUM OF THE SPOOL CP 
AND THE COLSYS CP. 


e 

NSUB 

CP 

TOTAL CP 

-2 

10 

84 

3.8 

3.8 

-4 

10 

168 

21 .3 

25.1 

-6 

10 

322 

40.8 

65.9 


TABLE 12. EXAMPLE 2 WITH Y = 2 AND Xq(0) = -4.29. NUMBER OF 

SUBINTERVALS (NSUB) AND CP TIMES TO SOLVE THE PROBLEM 
BY COLSYS WITH CONTINUATION IN e FROM E = 10“^. TOTAL 
CP IS THE ACCUMULATED TIME FOR THE SEQUENCE. 
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course to follow is perhaps to use asymptotic and numerical methods in 
tandem. For example, a rough numerical solution could be obtained for several 
values of e whxch could then be used to suggest the form of an asymptotxc 
solutxon. The asymptotxc approxxmatxon could then be used to refxne the 
numerxcal solutxon, and so on. It xs also possxble that sxngular perturbation 
theory could be used to construct special methods that are approprxate for 
specxfxc problems as e.g., in Flaherty and Mathon [9] and Ascher and Weiss [2, 
3, 4]. 
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